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Abstract 

The  object  of  this  project  is  aimed  at  development  and  application  of  a  new  3D  wave  propagation  and 
modeling  method  in  complex  heterogeneous  media  using  one-way  wave  approximation  to  study  many 
outstanding  problems  of  regional  wave  propagation  in  the  context  of  nuclear  test  monitoring.  The  great 
advantages  of  one-way  propagation  methods  are  the  fast  speed  of  computation,  often  by  several  orders 
of  magnitudes  faster  than  finite  difference  and  finite  element  methods,  and  the  huge  saving  in  internal 
memory.  The  first  year’s  effort  is  concentrated  on  developing  theory  and  techniques  for  fast  calculation 
of  reflected  or  backscattered  waves,  namely  methods  of  synthetic  seismograms  using  the  concept  of 
one-way  wave  propagation  and  single  interaction.  The  research  effort  supported  by  this  contract  has 
produced  the  following  results  regarding  the  backscattered  field  modeling.  Wu  and  Huang  (1995a,  b) 
has  introduced  a  wide-angle  modeling  method  for  acoustic  waves  using  the  De  Wolf  approximation  and 
phase-screen  propagator.  Xie  and  Wu  (1995)  introduced  a  complex  screen  method  for  the  calculation 
of  backscattered  elastic  waves  under  small-angle  approximation.  Wu  (1995)  extended  the  wide-angle 
method  to  the  case  of  elastic  waves,  and  derived  the  relation  between  the  wide-angle  method  and  the 
screen  approximation  for  the  acoustic  case.  Numerical  examples  showed  excellent  results  of  the  method. 
Comparison  of  angular  dependence  of  reflection  coefficients  calculated  from  synthetic  seismograms  by 
this  method  with  the  theoretical  curves  showed  good  agreement  between  theory  and  synthetics.  The 
wide-angle  version  of  the  method  can  even  model  the  critical  and  post-critical  reflections.  From  the 
comparison  of  synthetic  seismograms  made  by  this  fast  method  and  by  finite  difference  calculations  the 
validy  and  applicability  of  the  method  has  been  demonstrated.  Even  the  screen  approximation  method 
which  involves  a  small-angle  approximation  for  the  medium- wave  interaction,  can  produce  satisfactory 
synthetic  seismograms  for  many  practical  cases.  The  wide-angle  version  of  the  method,  though  more 
computationally  intensive  than  the  screen  method,  will  have  a  wide  range  of  applications.  The  next 
year’s  effort  will  be  concentrated  on  solving  the  free  surface  and  Moho  reflections  and  applying  the 
method  to  the  synthetic  seismograms  of  regional  phases. 
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1  Objective 


The  new  task  for  monitoring  a  Comprehensive  Test  Ban  Treaty  and  the  current  Nuclear  Non-Proliferation 
Treaty  presents  a  great  challenge  to  the  existing  methods  using  regional  phases.  This  is  because  the 
nature  and  properties  of  regional  phases  remain  unclear  due  to  the  lack  of  theory  and  efficient  methods 
for  computing  synthetic  seismograms  in  complex  environments.  The  object  of  this  project  is  aimed  at 
development  and  application  of  a  new  3D  wave  propagation  and  modeling  method  in  complex  heteroge¬ 
neous  media  using  one-way  wave  approximation  to  study  many  outstanding  problems  of  regional  wave 
propagation  in  the  context  of  nuclear  test  monitoring.  The  great  advantages  of  one-way  propagation 
methods  are  the  fast  speed  of  computation,  often  by  several  orders  of  magnitudes  faster  than  finite 
difference  and  finite  element  methods,  and  the  huge  saving  in  internal  memory.  The  new  project  will 
take  the  advantages  of  the  recent  progress  in  elastic  one-way  wave  propagation  theory  (Wu,  1994)  and 
continue  to  develop  the  theory  and  method  for  the  purpose  of  regional  wave  synthetic  seismograms. 
The  major  effort  for  the  new  project  will  be  to  solve  the  problems  of  backscattering  calculation  and 
free-surface/Moho  reflections  so  that  the  one-way  wave  theory  and  method  can  be  applied  to  the  re¬ 
gional  waveguide  environment.  Then  the  influence  of  upper  and  lower  crustal  small-scale  heterogeneities 
to  Lg  propagation  and  attenuation,  and  other  path  effects  will  be  examined  by  numerical  simulations 
using  the  new  method.  The  development  of  the  theoretical  model  will  be  combined  with  analysis  of  the 
Eurasian  regional  phase  observations  to  extract  quantitative  information  from  regional  phases. 

2  Research  Accomplished 

The  recent  successful  extension  and  applications  of  one-way  elastic  wave  propagation  methods,  e.g.  the 
complex  screen  method  (Wu,  1994;  Wu  and  Xie,  1994),  stimulated  the  research  interest  in  develop¬ 
ing  similar  theory  and  techniques  for  reflected  or  backscattered  wave  calculation,  namely  methods  of 
synthetic  seismograms  using  the  concept  of  one-way  wave  propagation  and  single  interaction.  The  re¬ 
search  effort  supported  by  this  contract  has  produced  the  following  results  regarding  the  backscattered 
field  modeling.  Wu  and  Huang  (1995a,  b)  has  introduced  a  wide-angle  modeling  method  for  acoustic 
waves  using  the  De  Wolf  approximation  and  phase-screen  propagator.  Xie  and  Wu  (1995)  introduced  a 
complex  screen  method  for  the  calculation  of  backscattered  elastic  waves  under  small-angle  approxima¬ 
tion.  Wu  (1995)  extended  the  wide-angle  method  to  the  case  of  elastic  waves,  and  derived  the  relation 
between  the  wide-angle  method  and  the  screen  approximation  for  the  acoustic  case.  The  validity  of 
the  method  and  the  wide-angle  capability  for  the  dual-domain  implementation  are  demonstrated  by 
numerical  examples. 

2.1  Synthetic  Seismograms  in  Heterogeneous  Media  By  One- Return  Approxima¬ 
tion 


A  new  method  based  on  multiple-forescattering  single-backscattering  (MFSB)  approximation,  i.e. 
the  one-return  approximation  (the  De  Wolf  approximation)  for  calculating  backscattered  fields  was 
introduced  (Wu  and  Huang,  1995a,  b;  Wu,  1995).  When  discontinuities  inside  a  medium  are  not 
very  sharp  or  parameter  perturbations  of  heterogeneities  are  not  very  strong,  reverberations  between 
heterogeneities  or  resonance  scattering  may  be  neglected.  However,  accumulated  effect  of  forward 
scattering  usually  can  not  be  neglected.  In  such  cases,  the  Born  approximation  is  not  valid  but  the  De 
Wolf  approximation  can  be  applied. 

The  Lipmann-Schwinger  equation  has  a  formal  solution  of  multiple  scattering  Born  series.  The 
widely  used  Born  approximation  is  the  leading  term  of  the  series.  After  renormalization  of  the  multiple 


532 


scattering  series,  De  Wolf  (1971,  1985)  derived  a  MFSB  approximation  given  by 

P(x)  =  pf(x)  +  k2  f  d3x'g^  (x;  x!)F(x!)pf  (x!)  ,  (1) 

J  V 

where  p!  and  g*  are  the  renormalized,  multiple  forescattered  field  and  Green’s  function,  respectively,  F 
is  the  perturbation  function  or  equivalent  force  which  is  a  scalar  function  in  the  case  of  scalar  media,  but 
involves  vector  operators,  such  as  gradient  and  divergence  operators,  and  tensor  operators,  such  as  the 
strain  operators  and  double  divergence  operator,  in  the  case  of  acoustic  and  elastic  media,  respectively. 
The  renormalized  pf  and  gf  will  be  calculated  using  the  phase-screen  or  complex-screen  propagator. 

To  speed  up  the  calculation  of  backscattered  fields,  the  local  Born  approximation  can  be  used  within 
a  thin-slab  (cf.  Fig.  1).  This  means  that  the  forescattered  field  pf  can  be  kept  unperturbed  and  g* 
can  be  replaced  by  a  constant  medium  Green’s  function  within  the  slab.  Assume  z'  and  z\  as  the  slab 
entrance  (top)  and  exit  (bottom)  respectively  (cf.  Fig.  1),  and  Fourier-transform  above  equation  with 
respect  to  xt,  we  get  a  dual-domain  expression  for  the  calculation  of  scattered  fields.  In  the  case  of 
scalar  media,  the  dual-domain  expression  is 

P(Kt,  z*)  =  ^k2  £  dze ^*‘-*1  J [  d2KTe-iK^(F(xT,  z)pf(xT,  z)]  .  (2) 

where  z*  is  the  receiver  depth,  xt  is  the  horizontal  coordinates,  Kx  is  the  horizontal  spatial  frequency, 
and  7  =  \Jk2  —  Klf.  Note  that  the  two  dimensional  inner  integral  is  a  2-D  Fourier  transform.  Therefore, 
the  dual-domain  technique  can  be  used  to  implement  the  equation. 

For  the  case  of  acoustic  media,  the  dual-domain  formulation  is 

P(KT,z*)  =  ^k2  J’1  dze^z*-z\  y  d2xTe~iKT'XT [eK(xT i  z)p* (xT,  z)} 

J  d2xTe~lKT'XT[£p(xT,z)Vpf(xT,z)\  J  ,  (3) 

where  eK(x)  =  ko/k(x)  —  1  ,  ep(x)  =  po/p(x)—  1,  and  k  =  |(Kt,  kz),  with  kz  =  ±7  for  forescattering  and 
backscattering,  respectively.  When  the  receiving  level  is  at  the  bottom  of  the  thin-slab  (forescattering), 
z*  =  zi; while  z*  =  z'  is  for  the  backscattered  field  at  the  top  of  the  thin-slab.  The  total  transmitted 
filed  at  the  slab  bottom  can  be  calculated  as  the  sum  of  the  forescattered  field  and  the  primary  field. 

The  dual-domain  expressions  for  the  case  of  elastic  media  is  more  complicated  and  involved  tensor 
operations.  Readers  are  referred  to  Wu  (1995). 

2.2  The  screen  approximation 

For  some  special  applications,  the  synthetics  only  involve  small-angle  backscattering.  In  this  case  the 
screen  approximation  can  be  applied  to  accelerate  the  computation. 

Under  small-angle  scattering  approximation,  we  can  compress  the  thin-slab  into  an  equivalent  screen 
and  therefore  change  the  3-D  spectrum  into  a  2-D  spectrum.  Dual-domain  implementation  of  the  screen 
approximation  will  make  the  modeling  of  backscattering  very  efficient.  In  the  following  the  case  of 
acoustic  media  will  be  shown  as  an  example.  For  the  screen  approximation  in  the  case  of  elastic  media, 
readers  are  referred  to  Xie  and  Wu  (1995)  and  the  forthcoming  publications. 

With  the  screen  approximation,  scattered  field  can  be  calculated  as 

P(Kt,  z*)  «  ij-e*  J [  dxTe~iKT'XT  S(xT)p°  (xT)  (4) 
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where 


rA  z 

5 (xr)  =  5v(xr)  =  /  dz[£K(x7',z)  -  £p(xr,z)]  /or  forescattering  (5) 

is  the  velocity  screen  and 

J/-Az 

d2el2fcz[eK(xT’,2)  +  £p(xt’,2)]  for  backscattering  (6) 

o 

is  the  impedance  screen. 

2.3  Numerical  examples 

First  we  compare  the  calculated  reflection  coefficients  by  the  one-return  method  with  the  theoretical 
ones  for  plan  wave  incidences.  The  model  is  shown  in  Fig.  2,  where  a  plane  interface  between  two  layers 
is  located  at  the  depth  of  500m.  The  model  is  defined  on  a  2048x300  rectangular  grid.  The  grid  spacing 
in  the  horizontal  direction  is  8m  and  that  in  the  vertical  direction  is  5m.  A  pressure  point  source  is 
applied  at  the  center  of  the  upper  border  of  the  model.  The  velocity  and  density  of  the  upper  layer 
are  2000m/s  and  1.0  g /cm3,  respectively.  The  frequency  range  used  in  the  calculation  is  from  14.6Hz 
to  19.5Hz  with  11  frequency  components.  We  calculate  the  reflection  coefficient  for  each  frequency 
component  and  take  an  average  over  the  11  frequencies. 

Fig.  3a  shows  the  results  for  10%  of  velocity  and  density  perturbations  in  the  lower  layer  and  Fig.  3b 
is  for  20%.  In  both  figures,  the  dashed  lines  are  theoretically  predicted  reflection  coefficients  of  plane 
wave  incidence  and  the  solid  curves  are  calculated  results.  We  can  see  the  good  agreement  between 
the  theoretical  ones  and  the  synthetics  when  the  incidence  angles  are  smaller  than  the  critical  angles. 
When  the  incidence  angles  are  near  and  beyond  the  critical  angles,  the  numerical  results  deviate  from 
the  theoretical  curves.  This  may  be  due  to  several  reasons  such  as  the  curved  wavefronts  of  waves  from 
the  point  source,  the  wavenumber  filtering  in  the  process  of  forward  propagation  which  reduces  the 
amplitudes  of  large-angle  scattered  waves,  and  the  effect  of  finite  layer  thickness.  When  the  velocity 
and  density  perturbations  of  the  lower  layer  are  -10%  and  -20%,  the  corresponding  results  are  given  in 
Fig.  4a  and  4b,  respectively.  We  see  from  these  figures  that  the  synthetic  reflection  coefficients  agree 
well  with  the  theoretical  results  when  the  incidence  angles  are  smaller  than  approximately  70°.  For 
larger  incidence  angles,  the  amplitudes  of  reflection  coefficients  decrease  because  of  the  wavenumber 
filtering. 

The  next  two  examples  is  for  the  test  of  screen  approximation  for  synthetic  seismograms.  Fig.  5 
shows  the  model  space  for  scalar  wave  modeling.  A  high  velocity  cylinder  with  10%  perturbation  is 
situated  at  the  middle  of  the  model  space.  A  plane  wave  is  incident  vertically  down  to  the  model 
space.  The  source  time  function  is  a  Gaussian  first  derivative  with  the  dominant  frequency  of  30 Hz. 
The  screen  interval  is  10  m  inside  the  cylinder;  while  the  horizontal  grid  spacing  is  3.125m.  Fig.  6a 
shows  the  comparison  of  synthetic  seismograms  calculated  by  the  screen  approximation  and  by  the  finite 
difference  method.  In  the  figure  only  the  transversal  components  are  shown.  We  see  good  agreement 
between  the  results  of  finite  difference  calculations  and  this  method. 

The  model  for  the  elastic  wave  scattering  is  similar  to  that  in  Fig.  5.  The  parameters  for  the 
background  media  are  a  =  3500 m/s,  (3  =  2050 m/s  and  p  =  2200 kg/mf.  The  inclusions  have  a  5% 
perturbation  for  both  P-  and  S-wave  velocities.  The  diameter  of  the  elastic  cylinder  is  300m.  A  plane 
P-wave  is  incident  on  the  cylinder.  Fig.  6b  gives  the  synthetic  seismograms  from  a  receiver  array  700m 
away  from  the  center  of  the  cylinder.  The  synthetics  marked  with  ’’SCREEN”  is  from  the  complex- 
screen  method  while  with  ”FD” ,  from  the  finite-difference  method.  The  upper  panel  is  for  x— component 
(transverse  component)  and  the  lower  panel  is  for  z— component  (longitudinal  component).  Shown  in 
the  figure  are  two  P-arrivals  and  two  S-arrivals  which  are  reflections  from  both  the  upper  and  lower 
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boundaries  of  the  cylinder.  The  results  show  general  agreements  in  both  amplitude  and  arrival  time. 
Since  the  receiver  array  is  only  550m  from  the  border  of  the  cylinder,  the  receiving  profile  spans  a 
rather  wide  range  of  scattering  angle.  That  means  that  even  the  small  angle  approximation  can  give 
satisfactory  results  in  many  cases. 

3  Conclusions  and  Recommendations 

We  have  made  significant  progress  in  modeling  backscattering.  Synthetic  seismograms  showed  good 
agreement  with  finite  difference  calculations  in  our  numerical  examples.  Next  year’s  effort  will  be 
concentrated  on  solving  the  problems  of  free  surface  and  Moho  reflections  and  applying  the  method  to 
simulate  regional  wave  propagations 
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Figure  1:  Geometry  of  the  thin-slab  for¬ 
mulation. 


•  x=(xT,  z) 


Pressure  Source 

.  \ 


Figure  2:  The  layered  model  used  in  the 
numerical  simulation. 


Reflection  Coefficient 


(a)  (b) 

Figure  3:  Comparison  between  theoretical  and  simulated  reflection  coefficients  for  a  high-velocity  ‘half-space’ 
shown  in  Figure  2  when  the  velocity  and  density  perturbations  of  the  lower  ‘half-space’  are  (a)  10%  and  (b)  20%. 
Dashed  lines  represent  the  theoretical  reflection  coefficients  and  the  solid  curves,  the  simulated  results. 


Figure  4:  Comparison  between  theoretical  and  simulated  reflection  coefficients  for  a  low-velocity  ‘half-space’ 
shown  in  Figure  2  when  the  velocity  and  density  perturbations  of  the  lower  ‘half-space’  are  (a)  -10%  and  (b) 
-20%.  Dashed  lines  represent  the  theoretical  reflection  coefficients  and  the  solid  curves,  the  simulated  results. 
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Phase  screen  vs.  FD 


Figure  6:  Comparison  between  the  results  of  the  screen  method  and  the  FD  method  for  scalar  waves  (a)  and 
elastic  waves  (b). 
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